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We solve the Bethe-Salpeter equation for a system of a heavy quark-antiquark pair interacting 
with a screened linear confining potential. First we show the spinless QFT model is inadequate 
and fail to describe even gross feature of the quarkonia spectrum. In order to get reliable descrip- 
tion the spine degrees of freedom has to be considered. Within the approximation employed we 
reasonably reproduce known radial excitation of vector charmonium. The BSE favors relatively 
large string breaking scale jU ~ 350Mey . Using free charm quark propagators we observe that 
J is the only charmonium left bellow naive quark-antiquark threshold 2mc, while the all ex- 
cited states are situated above this threshold. Within the numerical method we overcome obstacles 
related with threshold singularity and discuss the consequences of the use of free propagators for 
calculation of excited states above the threshold. 
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1. Introduction 

Excited meson spectroscopy is a keystone experimental output, which is essential for under- 
standing of quark-antiquark interaction. The dynamics of its constituents is dominantly driven by 
solely known strongly interacting quantum field theory- quantum chromodynamics. The explana- 
tion of confinement of constituents and all other colored objects is one of the great challenge of 
theory of strong interaction. The confining forces between quark antiquark inside mesons should 
lead a large degeneracy which emerges from the spectra of the angularly and radially excited res- 
onances. As reported in [|l|, ^, |3[ ^ such degeneracy is observed in pp annihilation by the Crystal 
Ball collaboration at LEAR in BERN ||5|]. Similar can be deduced from heavy quarkonia produc- 
tion in e+e^ anihilations. Recently, BABAR, Belle and BESS experiments continue collection of 
various meson experimental data. 

Heavy quarkonia, as relatively simple bound state systems made of the quark and antiquark of 
the same flavour, attract many theoretical and experimental attentions ||^. Following ideas of Ref. 
[^], confinement in heavy flavor hadron sector has typically been associated with a linearly rising 
potential between constituents The spin degeneracy observed in heavy quarkonia spectrum tell 
us that the main confining part of the interaction should be largely spin independent. A various 
lattice fits leaded to various predictions for the potential between quark-antiquark states. A well 



known is Cornell parameterization of the static Wilson loop derived potential M, 10], 



V{r) = -a/r+or (1.1) 

which appeared to be suited for description of the first few excited mesons. 

In the absence of dynamical quarks the nonrelativistic static potential shows up a linear asymp- 



totic of the interquark potential. Recently in Refs. [ |11| , |12| , |13[ ] it has been found that the meson 
spectroscopy is better described by "confining" potential which is bounded from above. The po- 
tential is flattened in large distance due to the string breaking and the screening mass effectively 
describes the inverse of string breaking range. To this point the exponential potential can be re- 
garded as a screened version of the linear potential in the region important for a few lowest states. 
The screening effect - the flatness of linearly rising potential arises form the string breaking sce- 
nario: including light quarks into the game then the creation of the quark-antiquark pairs, i.e. pions 
and other light mesons is energetically favorable. The screening is necessarily included to explain 
observed high radially excited heavy mesons- without any doubt, the spectrum deviates from the 
linear Regge trajectory. Recall, the deviations from the linear Regge trajectory is expected in the 



light meson sector as well [ 14 ] 



On the other side screening of Coulomb potential is more subtle matter, which could be related 



with soft gluon mass generation [16|] through the Yang-Mills Schwinger mechanism. Actually 
the dynamical gluon mass generation is suggested by finiteness of the lattice gluon propagator in 
the deep infrared, which fact is apparently pronounced in the recent lattice data [[l7|]. Both scales, 
the first given by the inverse of soft gluon mass and the second- the screening mass characterizing 
string breaking, are independent in principle, albeit one naturally expects their similar size ~ ^qcd- 
Using two models based on Bethe-Salpeter equation we study the effect of retardation in the 
slopes of radially excited vector mesons. In one of the model we completely ignore the spin degrees 
of freedom in order to see the effects of retardation and effect of off-shellness of the constituents 
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scalar propagators in fully Lorentz invariant BSE approach. The second model is based on quite 
standard assumptions- the (anti)quarks have spin 1/2 with the scalar confining part of the interac- 
tion and effective vectorial one gluon exchange. 

For bound states which lie above the naive quark-antiquark threshold, the BSE kernel be- 
comes singular, which causes that many usual numerical treatments fails and/or become impossi- 
ble in practice. To authors knowledge, theoretical explanation of plethora experimentally identified 
heavy quarkonia masses has been obtained only in relativistic quantum mechanic [^] or by solv- 



ing various phenomenological 3D reduction of BSE Ill9|, gO|]. The main motivation is to fill this 
missing gap and thus present our BSE study of above threshold excited bound states as a necessary 
step towards our better understanding of excited mesons in the nature. 

2. Model with spinless quarks 

In traditional quantum mechanical wisdom the spin interaction between quark and antiquark 
inside the mesons is necessary for relatively small hyperfine splitting and thus neglection of the 
quarks spin one should reproduce Regge trajectory for orbitally and radially excited mesons as 



well. Such picture has been confirmed in refs [21, 23] where spinless quarks has been considered. 



In the first model we just consider spinless quarks interaction through the relativistic invari- 
ant generalization of exponential and Coulomb(Yukawa) potential (see the Section next for the 
discussion of the nonrelativistic limit). To accomplish this we solve the following BSE 



T{p,p) = -ij ^v(k,p,p)d^\k,p)r{k,p) (2.1) 

where G is two body propagators 



gP1(A:,P) = D{k + P/2,ml)D{-k + P/2 
D{k,M^) ^ 



,ml) 



—m^ — He 

where is constituent quark mass and the interaction kernel is chosen such that 



Vsiq) 
and effective OGE 



V = V, + V, 
C 



with q = k — p and where /i = /x^ ~ (0.1 — 0.3)mc for simplicity (nic is the constituent charm 
quark mass). Recall here that BSE considered here (and in fact in the next section as well) thus 



represent generalization of nonrelativistic model considered [[TJ, [Tl|] where the potential (p^), 
complemented by coulombic term and constant positive shift o/jJ., provides quantum mechanical 
model describing experimentally known charmonium spectrum. 
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The BSE has been solved after performing Wick rotation and integration over the 3d space 
angles. It transforms the 4-dimensional BSE into the 2-dimensional one which has been solved by 
the method of the iterations by the similar method described in the recent paper [23]. The 2d BSE 
we have actually solved reads 



/OO /'CO 
dk4 / dy.T{k,P)G2% + K,\ ; 

k2 



(2.3) 



C 



Ak^ P 



In 



k^ + p^- 2k4P4 - 2|k[ IpI + /^^ 
k^ + p^- 2K4P4 + 2|k| |p + 



{kl-P^/4 + mlf + klP^ 



(2.4) 



Let us only mention here that the kernel of BSE becomes singular for the bound states heavier 
then sum of constituents, which requires certain numerical care when is solved numerically. 

For expected values of the parameters C,g the ground state appears near , but allays bellow 
the threshold 3GeV . This s-state corresponds with degenerate experimental partners of rjc and Jy 
mesons. We attempted to fit the distribution of other known s-states states, e.g. by adjusting C,g,}JL 
to get i//^(25'),i/a(35), etc... In order to ensure the numerical stability we also varied integration 
volume (in p-space) and the number of points as well. In order to get optimal density of integration 
points the upper boundary A is introduced in the momentum integrations. Let us stress, the wave 
function is exponentially decreasing and we always take A such that cut contribution is largely 
negligible ensuring that the main role of A is to adjust the density of integrations points and not to 
regulate the finite integrals. More precisely, for the variable k we map single integration interval 
(0, 1) into (0,A) by simple rescaling (and (—1,1) into (— A,A) for variable ^4). 

The results are shown in figures § for different volumes and its general features do not change 
considerably when changing the parameters in a range ±30%. While it is not difficult to fit two or 
three states we have failed to reasonably describe the experimentally known s — states as a whole, 
the obtained discrepancy is tedious since one allays gets the distribution of the states considerably 
dense then one can deduce from the known experiments. For instance, decreasing C does not lead 
large decrement of the density of the states , instead of, the numerical stability is lost, see Fig. ^. 
Increasing one can shift the intercept, but the rest of the energy levels are only slightly changed 
(also the numerical value is limited from the above by some critical value, above which the BSE 
spectrum become continuous). 



3. Model II, Vector BSE for mesons 

In principle the framework of Schwinger-Dyson and Bethe-Salpeter equations in QCD should 
offer unique fully Poincare invariant generalization of well known quantum mechanical picture of 
quarkonia, however a solution is always incomplete due to the truncation of the equations system. 
For instance how to truncate the set of equations in order to reproduce Wilsonian loop potential 
still remains unclear. Due to this fact we rather phenomenologically estimate what should be the 
form of the BSE kernels here. 
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Figure 1: Solution of BSE for C = 6Gey*, ^ 0.5 * (2;r)2Gey2,jU = 0.15Gey,m,. = 1.5GeV and various 

cutoff A. V"^^'^"!^ inrrp.a^iP.Q from thp Ipft to thp. ricrht 




Figure 2: Showing instability. The same as in fig. 1 but for C = 1 .6GeV 



In what follows we consider "confining" interaction kernel of the form 

which couples as a scalar in the quark-antiquark scattering kernel in the meson BSE. The second 
considered part of the interaction kernel has the Dirac decomposition identical with the one gluon 
exchange ~ T^pV^y^.p,. 
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For clarity we write down the BSE completely here 



S-\q + P/2)x{p,q)S-\q + P/2) = -i J -^,Y^X{k,P)YvG^'' {k - q) 

- ij j2^X{k,Pmk-q); 

(3.2) 

In words, scalar functions and are Poincare invariant generalization of quantum mechan- 
ical potentials. Here, clearly G^^ represents Gluon propagator in Feynman like gauge, where the 
effective soft gluon mass jXg has been introduced. Like in the previous scalar case, the double pole 
scalar interaction Vs leads to regular exponential potential in the position space. Actually, in heavy 
quark limit one can consider three dimensional potential 



Jo k 

-871 



(k2 + ^2)2 

where the potential in position space reads 

-fir 



(3.3) 



V{r) = -o- . (3.4) 



S in Eq. ( |3.2| ) stands for charm quark propagator, which in the simplest approximations reads 

S-\l)=f-mc, (3.5) 
and Xv represents the Bethe-Salpeter wave function which has the general form: 
Xv{q,P) =y^Zyo+ ^S.qXvi+ ^S.qXv2 + S.qXv3 

+ [^, /']ZV4 + [i^, Mv5 + f]XV6 + 175 fXvi , (3.6) 

with = e^va/3^''^"£'^, = -1, = 0. 

Munczek and Jain [ p5| ] shown that VO component is dominant for the all ground state mesons, 
which dominance is particular for mesons made form heavy flavor (anti)quarks. We assume the 
same apply for excited states as well and we neglect all other components in presented study. 

Within the approximation the Bethe-Salpeter equation in the rest frame reads 

XV0{PE,P) = , 2 2^M2IA\L ' (3-7) 



where 



h = -J ^[-2V,. + V,]xvoikE,P) ; 



{2n) 
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C 



(3.8) 



where we have performed Wick rotation into the EucUdean space, for Euclidean momenta kE = 
{k4,k), k\ = k\ + k^, while the total momenta is kept timelike = — P| = as required for 
the bound states. Again an extra problem arises for bound state which are heavier then sum of 
constituents quark masses. As we are using propagators with single real poles, the BS wave func- 
tion becomes singular since being proportional to the product of quark propagators. For the vertex 
function, the threshold like singularity should appear for the solution with > 4m^. Due to the 
presence of the kernel singularity more or less standard matrix methods [^] fail since the inversion 
of numerical matrices is not possible. Also we do not explore more or less conventional expansion 
into the orthogonal polynomials which loses its efficiency when, as one expects, relatively large 
number of polynoms is necessary. Instead of, we rather solve the full two dimensional integral 
equation by the method of simple iterations. For this purpose we discretize P^ and step by step we 
are looking for the solution of the BSE with given P^. Performing several hundred iterations for 
each PI we identify the solutions as those for which difference between iterations vanishes. 

Without loosing the generality the BSE can be transformed into the 2dim integral homoge- 
neous equations. It requires normalization. Instead of using physical one, to achieve a good numer- 
ical stability of the iteration process we implement normalization condition through an auxiliary 
function A(P). We choose A(P) such that it makes BSE nonlinear but mainly numerically stable 
and such that the BSE solution has been identified when A(P) = 1. Simultaneously the integrated 
difference a between consecutive iterations must vanish at the same time. We found that these two 
conditions happen simultaneously, while for other values of parameters P, X (P) 7^ 1 the numerics 
do not provide vanishing difference between iterations. For interested readers the details of the 
numerics is described in the paper [|23|]. 

As in the previous case of spinless BSE the numerical convergence of A — )• 1 and the resulting 
numerical solution itself crucially depend on the density of integration point. Therefore, similarly 
to the previous treatment we introduce cutoff A for the purpose. 

For clarity we restrict ourself to the dynamics of charmonia in this paper. The numerical values 
of the models are the following 

C = 5A\%GeV'^ ;a, =gV(4?r) = 0.2 (3.9) 
mc = \.56\5GeV;ix=lXg = 36A35MeV 

and we use the mass of meson to fine tune the correct scale at the end (m^ = 1 .5 ; C = 5 in units 
where M(7/*P) = 2975). The results presented in the Tab. |]have been calculated for A = lOOOGeV 
where the numerical stability has been (ahnost) achieved. We also compare with the experimental 
data in the Tab. [l] and Fig. ^. 

As one can see in the Table [T[ there aie two more BSE solutions with masses situated in 
between ^l" and i/a , while the rest of calculated excited states quite nicely agree with the data. We 
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Figure 3: Comparison of BSE solution with PDG data for vector charmonium. 



have not found possible way (by varying parameters C,a,mu,mc) to exclude these two additional 
states from the solutions. 

4. Discussion of the results 

We have formulated two BSEs models for vector quarkonia. Ignoring the spin degrees of 
freedom and treating quarks as spinless particles is not justified approximation when using BSE. 
While the presented numerical results do not represent ultimate numerical search we argue here 
that the retardation effects can be largely different in spinless and spin-incorporated BSE. In the 
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Table 1: Comparison of calculated BSE spectrum with PDG data (third column). The first column represents 
preliminary results presented at the talk at QCD-TNT2011. They have been calculated on the grid with 
A = lOOOGeV and 64*128 number of points. The second column represents the new results published in 
which has been obtained with 96*192 points. Quantum numbers correspond with assumed quantum 
mechanical assignment []l3[]. 
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second model -with two aforementioned exceptions of two additional states- the resulting spectrum 
is comparable with the experiments whenever the experimental data are available. The agreement 
between our results for higher states and the one measured in the experiments is impressive, re- 
maining difference between theory and experiments is due to the approximations, e.g. due to the 
interplay of quark and (D,D*) mesonic degrees of freedom, such couple channel effects are difficult 
to incorporate into the Bethe-Salpeter analysis presented here. 

At given stage the question of confinement is beyond scope of the presented work, however 
we expect some changes when confinement is correctly incorporated. First of all, we expect the 
quark propagators should not have a free particle pole and therefore the BSE could not posses 
ordinary threshold singularity in this case. In the paper we are dealing with BSE where the propa- 
gators describe free -instead of confined- quarks. Therefore the threshold singularity unavoidably 
appears as an artefact here. The lowest lying excited charmonia are the one closest to the naive 
quark-antiquark threshold and we naturally must expect some defects in the calculated spectrum. 
We argue two more states which appear are the artefact of inappropriate usage of the free quark 
propagators. Numerical solution and an implementation of the confinement at this level remains 
for future study. 

Furthermore, we expect that the knowledge of off-shell behavior (^4 dependence) of the BSE 
amplitudes can be important in various hadronic processes. Due to this it is worthwhile to study 
not only the static property like the spectrum here, but the cross sections of a processes including 
studied mesons. 
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